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Description 

Optimized calibration data based 
methods for parallel digital feedback, 
and digital automation controls 

Cross Reference To Related Applications 

[0001] This application claims the benefit of U.S. Provisional Ap- 
plication Ser. No. 60/408,92 lfiled Sep. 7 2002. 
Background of Invention 

Field of the Invention 

[0002] The present invention relates generally to design of multi- 
ple inputs multiple outputs (MIMO) system's control that 
allows nearly optimal operation over wide range of system 

states. 
Background 

[0003] Control design techniques used for achieving nearly opti- 
mal control solutions (1) are based on knowledge of inter- 
nals of the system (further referred as a plant). Such 
knowledge is acquired from analytical models of a plant or 



training process that assumes no knowledge of the plant's 
internals (2,3)- Both approaches are well known to engi- 
neering community and both represent significant mile- 
stones on a way of design and implementation of nearly 
optimal control solution. 
[0004] Analytical approach represents less challenge since it is 
based on explicit knowledge of a plant and its possible 
states. Art of design of control solution in this scenario 
narrows down to solving model equations to select robust 
and at the same time efficient solution. Robustness and 
efficiency contradict each other. Robustness addresses 
stability of a control solution at variations of model equa- 
tions, while reducing efficiency of control that make it less 
optimal. 

[0005] Training approach experience quick evolution and is in 
focus of interest of many scientific and engineering 
schools. This approach assumes no prior knowledge of a 
plant and plant model is discovered using real-time mon- 
itoring of a plant control variables at various control in- 
puts and auxiliary parameters values (4). This approach 
requires little or no skills from design engineer and 
tempts to discover allowable domain of a plant states and 
plant transfer functions. Complexity of this approach is 



bound to algorithm of such discoveries. Some of such al- 
gorithms use fuzzy logic and neural networks to build in- 
ternal state model of a plant. Main factor in selection of 
training approach is speed and stability of the model. 
While allowing model to be easily mutable makes control 
solution less robust with respect to unexpected distur- 
bances, making solution to too conservative increases 
length of training process and could skip discovery of 
volatile system states. Control design engineer's choice is 
also driven by considerations of solution stability in case 
of partial plant failure (5). This factor usually rules out an- 
alytical approach as incapable of adaptation to significant 
model changes. Contrary highly mutable modeling algo- 
rithms address such events with highest efficiency. 
[0006] The purpose of the training is to discover a plant's im- 
pulse transfer functions. This task usually solved through 
monitoring of the plant's response on a control input vec- 
tor. Often these vectors are selected to have shape of 
Heaviside step function. Downside of this approach 
emerges from uneven power spectrum of the transfer 
function. Measurement of noise power spectrum reveals 
significant variations of S/N ration across desired fre- 
quency range. Ideal solution for this problem would be 



selection of input vector with power spectrum that com- 
pensate for described variations. Causality principle pre- 
vents such selection due to unknown plant's dynamics. In 
current state of the art this causes the problem of non 
uniform model accuracy across operational frequency 
range. Current Invention addresses this problem. 

[0007] Another factor that only recently became possible to 

model analytically is hysteresis. Hysteresis imposes limi- 
tations on majority of LTI based models that either ob- 
tained analytically or result of real-time plant discovery. 
This phenomenon also increases complexity of neural 
network models (6). 

[0008] Another general aspect of current state of the art that re- 
lates to current invention is selection of time-optimal 
control functions for MIMO plants. This subject was ad- 
dressed in great details by Timothy D. Tuttle and Warren 
P. Seering (1). They propose elegant solution for discovery 
of time-optimal commands for linear systems. Neverthe- 
less suggested solution for this problem remains hard to 
implement in cases of saturation of controller outputs and 

it is hardly applicable to cases of dynamic plant's state. 
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Brief Description of Drawings 

[0016] The present invention will be more readily understood 

from a detailed description of the preferred embodiments 
taken in conjunction with following figures. 

[0017] FIG. 1 is qualitative comparison of spectral composition of 
a control response and random disturbances and mea- 
surement noise. 

[0018] FIG2. is a recursive algorithm of system identification with 
multiple steady states. 

[0019] FIG. 3 is a schematic use of the invention as add-on com- 
ponent responsible for validation and refinement of 
plant's response model. 

[0020] FIG4. is an example of random and systematic approach 
to discovery of plant's hidden states. 

[0021] FIG. 5 is an example of constraint's space segmentation. 

[0022] FIG. 6 is a simplified diagram of two channel MIMO con- 
troller. 

[0023] FIG. 7 is an illustration of optimal transition trajectory for 
plant's output. 

[0024] FIG. 8 is illustration of close-loop feedback operation. 

[0025] FIG. 9 is a block diagram of multithreaded MIMO con- 
troller. 



[0026] FIG. 10 is a block diagram of calibration channel. 
[0027] FIG. 11 is a block diagram of feedback channel. 
[0028] FIG. 12 is a block diagram of feed -forward channel. 
Detailed Description 

Common Conversions 

[0029] a method and apparatus for optimal digital control of var- 
ious processes and its use as a close-loop feedback and 
feed-forward system is described. In the following de- 
scription, numerous implementation details are omitted. It 
will be apparent, however, to one skilled in the art, that 
the present invention may be practiced without those spe- 
cific details. In most instances, well-known structures and 
devices are shown in block diagram form, rather than in 
detail, in order to avoid obscuring the present invention. 

[0030] Some portions of the detailed descriptions that follow are 
presented in terms of algorithms and symbolic represen- 
tations of operations on data bits within a computer 
memory. These algorithmic descriptions and representa- 
tions are the means used by those skilled in the data pro- 
cessing arts to most effectively convey the substance of 
their work to others skilled in the art. An algorithm is 



here, and generally, conceived to be a self-consistent se- 
quence of steps leading to a desired result. The steps are 
those requiring physical manipulations of physical quanti- 
ties. Usually, though not necessarily, these quantities take 
the form of electrical, magnetic signals, location, velocity, 
acceleration of physical body, temperature or thermal 
waves capable of being measured, stored, transferred, 
combined, compared, and otherwise manipulated. It has 
proven convenient at times, principally for reasons of 
common usage, to refer to these signals as bits, values, 
elements, symbols, characters, terms, numbers, or the 
like. On other hand all these quantities referred as re- 
sponses, results, outputs, or the like when they consid- 
ered as a physical characteristics that are subject of con- 
trol. 

[0031] it should be borne in mind, however, that all of these and 
similar terms are to be associated with the appropriate 
physical quantities and are merely convenient labels ap- 
plied to these quantities. Unless specifically stated other- 
wise as apparent from the following discussions, it is ap- 
preciated that throughout the present invention discus- 
sions utilizing terms such as "processing" or "computing" 
or "calculating" or "determining" or "displaying" or the 



like, refer to the action and processes of a computer sys- 
tem, or similar electronic computing device, that manipu- 
lates and transforms data represented as physical 
(electronic) quantities within the computer system regis- 
ters and memories into other data similarly represented as 
physical quantities within the computer system memories 
or registers or other such information storage, transmis- 
sion or display devices. 
[0032] The present invention also relates to apparatus for per- 
forming the operations herein. This apparatus may be 
specially constructed for the required purposes, or it may 
comprise a general-purpose computer selectively acti- 
vated or reconfigured by a computer program stored in 
the computer. Such a computer program may be stored in 
a computer readable storage medium, such as, but is not 
limited to, any type of disk including floppy disks, optical 
disks, CD-ROMs, and magneto-optical disks, read-only 
memories (ROMs), random access memories (RAMs), 
EPROMs, EEPROMs, magnetic or optical cards, or any type 
of media suitable for storing electronic instructions, and 
each coupled to a digital processing device. The algo- 
rithms and displays presented herein are not inherently 
related to any particular computer or other apparatus. 



Various general-purpose machines may be used with pro- 
grams in accordance with the teachings herein, or it may 
prove convenient to construct more specialized apparatus 
to perform the required method steps. The required 
structure for a variety of these machines will appear from 
the description below. In addition, the present invention is 
not described with reference to any particular program- 
ming language. It will be appreciated that a variety of pro- 
gramming languages may be used to implement the 
teachings of the invention as described herein. 
[0033] Some portions of the detailed descriptions use terms 
transfer and domain to represent known mathematical 
operator that convert entities that represent some physi- 
cal quantities into different forms. Such transformations 
do not alter meaning of physical entities while transferring 
them to different representation domain. It should be 
borne in mine that it is obvious to one skilled in mathe- 
matics and signal processing that there is equivalence be- 
tween representations of physical entities in different do- 
mains. Keeping this in mind to prevent repetitive descrip- 
tions of equivalent procedures performed using alterna- 
tive transforms following detailed description uses Fourier 
transform only, which in most cases can be replaced by 



wavelet transform or Z transform or Laplace transform. 
Preferred Embodiments 

[0034] This invention discloses innovative approach to discovery 
of impulse response function of MIMO plant with linear 
time-invariant state's model. This approach benefits 
higher model accuracy across complete frequency spec- 
trum. Considering plant stimulus vector s(t) that exist in 
space s=Un Y, where uc% n is space of control inputs and 
Yc% m is space of plant's detectable outputs. The s(t) vec- 
tor of functions consists of commonly used vector of in- 
puts u(t) for the plant and vector of plants outputs y(t). In 
assumption of LTI plant [/and Y spaces are bound by 
plant response matrix G , so y(t)=g (t)*u(t)(eq.l0l), where 
matrix g (t)G% mxn known as impulse response of the 
plant. This equation can be extended to plants with dy- 
namic state models only with restriction that state of the 
plant remains steady for complete subspace of u(t) and 
y(t). 

[0035] jo simplify description of the invention it is assumed that 
in its current state plant behave as LTI system when space 
of u(t) and y(t) restricted by small sphere O, where Oeft m+n 
and has its center at (u(0),y(0)). In such case inverse rela- 
tionship between u(t) and y(t) can be solved u(t)=g (t)* 



y(t)(eq.l02), which represents feed forward term of a con- 
troller. Matrix g i (t)G% nxm . 

[0036] p| an t input vector u(t) consists of n scalar functions. Cur- 
rent control theory considers these functions as indepen- 
dent. This assumption is only a special case of more 
generic approach that considers cross-correlations be- 
tween different inputs. Deconvolution matrix g ft) for 
vector of plant inputs can be found through well-known 
computational and analytical techniques, so cross depen- 
dency of inputs could be defined asu(t)=g (t)*u(t)(eq.l03). 

[0037] LTI restrictions impose limitation on plant output vector 
y(t). Cross-dependency of individual outputs defines ma- 
trix g (t) that remains invariant for any value of s(t)eO 
vector. These relationships result in equation y(t)=g (t)* 
y(t)(eq.l04),\N\\\ch describes plant behavior in absence of 
control inputs. This set of equation describes extended 
model of the plant and in combination with eq.103 de- 
scribes a system composed of a controller and a plant as a 
single entity, for convenience it will be referred as an ex- 
tended plant. Vector r(t)e% m n represents vector of virtual 
output of an extended plant that is describe by r(t)=g(t)* 
s(t)(eq.l05), where s(t)=(u(t),y(t)) T and 
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[0040] 



k, ^- , (m+n)x(m+n) . , . 

. Notice that g(t)G% is square matrix. 

Example of power spectrum of a plant's impulse response 
compared with power spectrum of measurement noise 
and disturbances shown on Fig. 1. It is often the case to 
have Response power spectrum attenuating at higher fre- 
quency faster than disturbance spectrum. Measurements 
noise usually has much wide spectrum than plant's re- 
sponses to control inputs and disturbances. Such behav- 
iors cause significant variations in signal to noise ration 
for different bands during collection of the plant response 
data, which results in less adequate control model. 
This embodiment discloses method for reconstruction of a 
plant response that provides predetermined accuracy □((*>) 
across full range of frequencies 102. Block diagram of this 
method shown on Fig. 2. 

1) Stimulus vector for a plant or extended plant is selected 
based on causality principle so only independent compo- 
nent of the vector receive control inputs. Initial inputs s 
(t) are selected as scalar functions with wideband spec- 

0i K 



trum 101. This selection can be determined based on 
knowledge of the plant internals or generic set of func- 
tions like Heaviside or Ramp or Rectangular can be used 
as zero approximation. 

[0041] 2) Response vector r(t) is acquired for predefined time pe- 
riod T q and its power spectrum is computed 104. The 
time period for data collection can be determined based 
on knowledge of the plant dynamics, be dynamically de- 
termined using incremental transforms to frequency do- 
main from convergence of static gain vector, or using 
other approaches known to one experienced in the art. 
Convergence of zero frequency gains can be estimated 
using known statistical approaches. For simplicity Fig. 2 
shows the use of Kalman filter 105, it is obvious to one 
experienced in the art that other well known techniques 
can also be used without limiting the invention. 

[0042] increasing T q 109 increases accuracy of static gain detec- 
tion. It also provides accurate estimate for error covari- 
ance of model data at low frequency band. Acquisition 
process is stopped when error of static gain drops below 
□(0) 106. 

[0043] 3) |_ ow Danc | model analysis 107 finds [u> ~,u> + ], where 

k k 

error covariance of the model across different time win- 



dows exceeds limits defined by □(u>) 102. Set of new s 
(U) ~,u> + )(t) is selected for each spectral band. Addi- 

k k 

tional steps are taken to optimize this set and replace 
each pair of {s (to + )(t)s i (co ',uj + M 

k k k k+1 k+1 k+1 

when total length of execution of s and s exceeds ex- 

a k k+1 

ecution time of a single vector s(uj ,00 + )(t) 108. 

3 k k+1 

[0044] 4) Selected vectors {s } are applied to the plant in order of 

k 

increasing oo " 103, each vector repeats until the plant 

k 

model converges at oo "to satisfy d((jo ~) criteria. During 

k k 

this process plant model is revalidated and invalidated [go 

p 

~,(jo + ] are excluded when d([(jo ,oo + ]) criteria is met. 
p p p 

Additional steps are taken to satisfy LTI assumptions for 
stable plant state. These steps include cluster detection 
methods 110 for validated segments of the plant model. 
To keep description simple no discussion on available 
clustering algorithm will be made. It is obvious to one ex- 
perienced in the art that multiple algorithms are available 
and for the goal of this invention the details of those al- 
gorithms are not important. As an example simple Stu- 
dent T test can be used to detect clustering by single re- 
sponse component. Method continues on this step to next 
invalidated spectral band until all bands are validated or 
clustering is detected. 



[0045] 5) Detection of clustering indicated transition of the plant 
to a new distinct state g . Such transition contradicts ini- 

k 

tial assumption on state's stability. When such transition 
occurs an assumption is made that newly acquired state is 
also stable for new O G% m+n . Previously constructed 

k 

plant model g then stored 111 in completed or incom- 

k-l 

plete state and new model construction begins 112 fol- 
lowing the described steps. 
[0046] The advantage of method described in current embodi- 
ment is its independence from the plant state as well as 
control inputs. That makes it possible to apply current in- 
vention to an extended plant which model needs to be 
revalidated or tuned. The fact that all stimuli s(t) used in 
this method are limited by OeR m+n allows the method to 
be applied during operation of the plant. Notice also that 
all operations of the method are linear, so it can be ap- 
plied on a real physical plant or a virtual plant constructed 
from modeling errors of an existing extended plant 201. 
In this case considering the extended plant with periodic 
or repetitive states transitions, the method of invention 
can be modified as shown on Fig. 3 to allow use of control 
inputs of the plant as a vector of stimuli to perform recur- 
sive converging of distinct states models while the plant 



performs transitions between those states. Referring to 
Fig. 3 an extended plant 201 includes basic controller 
202. Add-on controller 203 operates according to the 
method of this invention. 
[0047] Although no robust method exists that is capable of dis- 
tinction between different steady states of a complex 
plant this subject remains in constant focus of the art. 
This invention proposes yet another method for solving 
this problem. It makes an assumption that no two states 
share the same stimuli vector s(t). Main restriction of this 
assumption is an existence of states transitions. It is al- 
ways possible to select an extended plant in such a way 
that the assumption will fail. Nevertheless for each ex- 
tended plant that is controllable or observable the as- 
sumption is valid. As an example a plant containing a mix 
of water and ice will behave as a plant with three distinct 
stable states (water only, ice only and the mix). Each state 
behaves in LTI fashion, but transitions between them may 
not be traced. If this is the case, consider temperature as 
an output of this plant. The temperature of the mix state 
could be equal to temperature of any other state, so the 
assumption mentioned above fails. Addition of pressure 
or volume or integral of supplied heat as second compo- 



nent to the response vector resolves the problem as the 
plant becomes observable. An addition of integral of sup- 
plied heat as a component of stimuli vector will make the 
plant controllable and also will make aforementioned as- 
sumption valid. 

[0048] Considering that an extended plant is observable or con- 
trollable, current stable state of the plant can be deter- 
mined by matching current s(t) to one of existing stable 
state model g using the eq. 105, where r(t)-s(t)eO .Such 

k k 

detection can be performed a priori or posteriori. A priori 
state transition performed as a part of dynamic feed for- 
ward control, while posteriori in dynamic feedback that may 
adjust feedback parameters based on current state of the 
plant. 

[0049] yet another embodiment of this invention discloses the 
use of special control inputs optimized for plants with 
hysteresis. The nature of a hysteresis described as alter- 
nation in plant model behavior due to changes in trajec- 
tory directions of state or input variables. This invention 
proposes the use of control inputs that can be repre- 
sented as finite functions with all or first k derivatives 
having the same sign over the length of the input. To per- 
form plant identification using such control inputs various 



well known methods of system identification can be used 
with some alterations. The use of method described in 
previous embodiment will be shown as an example. 
[0050] Considering control input having n components and each 
component having u + .(t) or u~ .(t) shape. There are 2*n 
input vectors required to perform complete hysteresis 
identification for each steady state. Advantage of using 
this approach is simultaneous detection of a hysteresis 
that could be caused by different derivatives of the input. 
Each of vectors u(t) has u(t) composed of components with 
opposite signs. The method described in the first embodi- 
ment can operate without changes by using random se- 
quence of proposed inputs. As an example consider tem- 
plate function 

f (/) = (e* - e~ l )j(t + 0.5) + H(0 

that can be scaled by time to satisfy frequency domain re- 
quirements of the method, and can be scaled and offset 
by magnitude to meet actuators constraints. In this case 
input vector u(t)=^E,(at+b)B+Q where B is signed vector of 
amplitudes and C is vector of offsets both selected to sat- 
isfy u(t)cO(s). 

[0051] The ultimate goal of the system identification methods is 



to uncover all hidden states of a plant in shortest time. 
This task can be accomplished using different strategies, 
some of which are: random probing of a plant state space, 
grid, etc. Fig. 4 illustrates these concepts. Random prob- 
ing 301 arbitrarily select B and C vectors to cover all sub- 
space of allowable actuator values and than builds map 
between the plant states and input vectors. Grid probing 
302 covers all space of inputs with uniform grid and per- 
form sequential or random travel through its nodes. Any 
selected strategy can use amplitude scaling to probe input 
space with higher density in regions adjacent to detected 
state transitions. 
[0052] plurality of LTI like models are constructed as a result of 
described processes and said models compose response 
model of the plant. This model contains explicit knowl- 
edge of the plant response at various states and various 
vectors of derivatives of the plant inputs. In absence of 
hysteresis at some state of the plant there will be only one 
LTI like model regardless of direction of vectors of deriva- 
tives of the plant inputs. In presence of hysteresis the 
same state will have multiple LTI like models where each 
model is associated with at least one direction of said 
vector of derivatives of the plant input. 



[0053] The application of this explicit hysteresis identification 
may include but is not limited to: generation of optimal 
control input solution for the plant that provide optimal 
trajectory of plant output transition; reduction of plant 
model that excludes hysteretic behavior by enforcement 
of additional constraints on control inputs. 

[0054] yet another embodiment of present invention discloses 
nearly time-optimal method of feedback control that 
comprises feed forward elements. Data of MIMO dynamic 
model of a plant enables one experienced in the art to im- 
plement methods for automated discovery of time-op- 
timal or nearly time-optimal inputs for transition of the 
plant from one steady state to another. For clarity of cur- 
rent disclosure no attempt to describe any of such meth- 
ods will be made since they are well described in refer- 
ences. It is important to notice that computation of such 
time-optimal commands for MIMO systems has high cost 
from digital resources perspective (processor utilization 
and memory use). Another important aspect to remember 
is dependency of such input model from various con- 
straints of an extended plant (such as actuators limits, 
power, fuel rate, etc.). These two factors make it unpracti- 
cal to implement true time-optimal commands in MIMO 



systems. Current embodiment describes method that 
keeps input commands close to time-optimal while per- 
mitting real-time implementation on MIMO plants that not 
necessarily confined to LTI. 
[0055] Constraints of an extended plant can be represented as a 

h 

vector of functions c (t,s(t))GC(t)c% for each steady state 

k 

of a plant. Such form of constrains indicate that they may 
change with time so volume of space C is dynamic but in- 
dependent of s(t). This embodiment assumes that rate of 
C(t) changes is small in comparison with rate of changes 
in s(t). If such assumption can not be made for some ex- 
tended plant then it is still possible to use this invention 
by selecting C =min(CA where operator min() takes region 
of space Cthat remains steady over characteristic period 
of s(t) changes. Because of imposed restriction it is safe to 
approximate C(t) with series of steps as 

c(o=2>,#(^) 

■ 

I 

, where D . is incremental change and H is Heaviside step 
function, and consider c(t) as a set of constant spaces C . 

i 

For each c . vector of constraints can be represented as c 



(s(t)). Important note here is that dimension of constraint 
vector may change with transitions between steady states 
of an extended plant. 
[0056] pig. 5 illustrates an example of constraint space C . For 

i 

simplicity of illustration the space boundaries 401 drawn 
as a rectangle and number of dimensions equals two. In 
reality boundary has more sophisticated multidimensional 
shape. This region of space is divided by plurality of 
nodes. The nodes can be selected using various desirable 
strategies. As an example on Fig. 5 the nodes are selected 
using grid method with linear dependency of the grid size 
from distance to the boundary and limited maximum and 
minimum size 




d = max(min( / , d m ), d^ n ) 



, where d and d are maximum and minimum grid 

max min 

sizes and r is shortest vector between current node and 
the boundary. It is obvious to one experienced in the art 
that choice of particular strategy affects number of nodes, 
optimality of control solution and other aspects such as 
allocation of digital resources, but it does not limits the 
invention. That is why no other examples of node selec- 



tion strategies will be made here. It is also obvious that 
different strategies can be used for various steady states 
of an extended plant. 

[0057] set of constraints for each node is well defined so time 
optimal control inputs can be computed and stored as 
time series of s values. These recipes are organized in 
lookup table that indexed by each dimension of the con- 
straint vector. Such indexing schema allows fast real-time 
access to recipes data. Operation of an extended plant 
encounter necessities of time optimal control inputs use. 
Each such event defined by current location of the plant in 
constraint space, which is known by history of previously 
issued inputs recipes and current s vector. Second defin- 
ing factor is static vector of desired state transition 
(switching distance), which represents only difference be- 
tween current and desired s vectors. 

[0058] Control input recipe then retrieved using indexed lookup 
to nearest node that closer to C boundary. Recipe con- 

i 

tains a set of solutions. Each solution is defined by maxi- 
mum switching distance and switching time. This set is 
ordered or indexed by maximum switching distance. So- 
lution with smallest switching distance that exceeds de- 
sired state transition is selected and scaled down to match 



the transition. Advantage of this invention is that it allows 
MIMO plant to be controlled during state transitions with 
desired control accuracy while remaining nearly time- 
optimal. 

[0059] Example of response spectrum reconstruction 

[0060] Following example illustrates the use of method of inven- 
tion to obtain desired model accuracy across required fre- 
quency range. For purpose of illustration it is assumed 
that matrix g and g are zero matrixes and space S sep- 
arated on two independent subspaces t/and Y. It is also 
assumed that both spaces have equal dimensions, and 
plant is LTI. 

[0061] | n accordance with all aforementioned conventions the In- 
vented method description considers equal number of the 
outputs and the inputs. In general case application or any 
changes in one of the inputs 



u 0 (t - 1 ) 



u p (t-t i )_ 

will cause changes in some or all outputs of the plant. A 

set of u (t) functions represents a set of stimuli simulta- 
p 

neously applied to the plant at the moment t . This pro- 

i 

cess provides valuable information about properties of the 
plant. The changes in the current state will continue to 
occur for some period of time after changes in the stimu- 
lus have stopped. These changes are acquired from the 
moment t . of initial change in the input and stored as a 
vector of response functions 



y,(f-0 



[0062] This data will be acquired for all inputs. Collection of vec- 
tors u forms a matrix of stimulus functions: 



S(0 = 



u° p (t-t 0 ) 



u l 0 (t-tj) 



u l p (t-t t ) 



[0063] Lower index p identified the stimulus component and up- 
per index / identifies a sequence the stimuli have been 
applied in. Collection of the response vectors forms the 
matrix of response functions 



yt(t-t 0 ) 



R(0 = 



yt(t-t 0 ) 



[0064] Lower index k identified the response component of the 
plant and upper index / identifies a sequence of recording 
identical to sequence of stimuli application. In order to 
complete calibration a matrix G of transfer functions of 
the plant is calculated from the equation: 

R = G*S 



R(t) = 



g°o(t) 



gl(t) 



yk ~ 2-1 

p 



* 



s(/) 



(eq.201) 

[0065] This equation is analogous to eq.101 with only difference 
that all components of input and output vectors are 
shown in group form. In the method of Invention it is as- 
sumed that the plant repeatedly passes through the same 
steady state. This example uses step function as a tem- 
plate for components of input vector: 



r o 



S(t) = 



QqH{1 — Iq) ... a^Hit — t^) ... 



alH(t-t 0 )...a l v H(t-t l ) 
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, where a is scalar coefficients and det(S(oo))=t=0. In this case 
the equation (eq.201) transforms into: 



F(R(0) = F(G(t))F 
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a 0 V ,v ' 0 



^0 

iG)F(R(t)) = F(G(0)A 



a\e mi 



a 1 e~ iwtl 
p 



■J/ 



\ 



, where co is a frequency 5 is delta function, H is heaviside 
step function. 

[0066] From this equation the matrix G of transfer functions can 
be obtained as 



F\icDF(R(t))A [ ) = G(t) 



-i 



(eq.202). 



[0067] | n case 0 f non square matrix generalized inverse can be 
used. In many real-world systems it is often the case that 
plant response has significant attenuation at some fre- 
quency range and their actual value become comparable 
with external interference and/or noise of measurement 
equipment. In this case results obtained from equation 
(eq.202) may become inaccurate. Use of repetitive or peri- 
odic input to improve model accuracy was disclosed in the 
method of Invention. This example demonstrates case pe- 
riodic input functions 

S («) = S(, - n ) 

(eq.203), where (jo q is period frequency. These stimuli ap- 
plied to the equation (eq.202) eliminates all components 
with lower frequencies and increase SNR (signal to noise 
ratio between response R and noise N-fold where N is 
number of full periods of S(t) executed during the calibra- 
tion. 

[0068] | n this case the following equation gives true improved 
G(t) matrix of transfer functions: 



, where u> is a frequency, w q =0, N h is number of whole 
periods used for calibration and B is Boxcar function. Re- 
constructed of complete response matrix 



G(/) 



z 



, where 



G;(o=F-fF(G t «))A t ^ ( y^ ) > 

(eq.204) 

[0069] This equation uses data received from static gain calibra- 
tion G o (u>) as well as data from calibration with periodic 
stimuli (eq.203) G (oo). Equation (eq.204) demonstrates 

h 

the case of calibration that uses any given number of 
stimuli with distinct periods. A for h^O is defined as: 

h 




(t-*l) & h 



In — , where IT is rectangular function 
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(eq.205) 

[0070] shown example illustrated benefits of plant model incre- 
mental reconstruction when accurate responses can not 
be obtained using traditional techniques. Example used 
rectangular functions as a stimuli and spectrum merge 
was done using boxcar functions. This choice was made 
for simplicity of illustration only and does not restrict the 
invention to choice of particular shape of the function. 

[0071 ] Example ofMIMO feedback controller with optimal control input 

[0072] Example illustrates operation of close-loop feedback con- 
troller that implements optimized set of control solution 
to achieve performance improvement. To facilitate subject 
illustration only necessary details of implementation are 
shown. 

[0073] The apparatus referred in this example schematically 



shown on diagram Fig. 6, for clarity of the demonstration 
it will have only two drivers that apply control input to the 
plant 601 and two sensors that characterize plant output 
or its state. Plant state can be changed so sensor A read- 
ings will follow curve R (t) and readings of sensor B re- 
main constant. To produce this type of results drivers 1 
and 2 should output D ^ (t) and D (t). 
[0074] Equation (eq.101) gives exact formula for computing 
these driver functions: 

F(rl) = ^F(g' k )F( S ' i ) 
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where symbol A indicates frequency domain 

R = GS => S = G-'R o F(sl) = £F(g*)F(ift 



where symbol represents an element of generalized inverse matrix 

z 
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(eq.206) 

[0075] Matrix G _1 can be obtained as a result of calibration pro- 
cedure as described in previous sections. 



[0076] 



A = F~' {F{g? W(R A )) + F 1 (F(g f )F(^ )) 



F\F(gf)F(R A )) + R B g 



D 2 =F l {F{gt )F(R A )) + F 1 {F{g B 2 W(R B )) 



R B =const 



F\F{g A 2 )F{R A )) + R B g 



In common case of N channels : 



D^t) = F-^r h ) + ^R k H 




(t)=const, which shows practical approach of eliminating 
coupling between different plant outputs. 
[0078] Matrix G -1 in frequency domain can be stored as a fre- 
quency series and partitioned by row. Each row represents 
separate stimulus/driver channel and independent of 
other channels. In general case of N channels equations 
(eq.206) defines pre-computed time series for each stim- 
ulus/driver. These time series can be partitioned by a 
driver channel and a set of desired response functions. 
[0079] such apparatus can efficiently operate in close-loop mode 
when readings from sensors contain disturbances from 
expected system state. Following demonstration shows a 
disturbance correction on cannel A. This is done just as an 
example and does not post limitation on the invention, 
the same approach equally applicable to any number of 
channels. At initial moment deviation of sensor A from 
nominal is -E Obviously, it is impossible to compensate 



a . 



this deviation momentarily because it will violate causality 
principle in events propagation. To compensate the devia- 
tion some correction curve should be selected that has an 
ending value equal E . It is reasonable to select a smooth 
trajectory to eliminate undesirable disturbances as shown 
on Fig. 7 with trajectory amplitude constrained by E . Al- 
though multiple approaches for selection of optimal and 
nearly optimal trajectories were described in listed refer- 
ences for simplicity of this example a curve defined as ViE 
erfc(-t/b) is considered, where 2b is effective halve width of 
the curve. In this case equations (eq.207) transforms to: 



F 



erf 



nk 



exp 



{-{nkbf): 



F 



( ( t \\ 
erfc -- 

. V bjj 
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A = E A F~ 



D 2 = E A F- 



where Ijzk = co 



(eq.208) 

[0080] N ote: erf(2)=0.9953, replacing erf(x) for |x|>2 with 1 and 
1 accordingly will introduce 0.23% relative error. Such er- 
ror is usually acceptable in majority of real-world applica- 
tions. These equations also demonstrate that decrease in 
response time of feedback causes an increase in power 



utilization of drivers. Because of physical limitation the 
driving power is limited. The equation (eq.208) allows ef- 
ficiently calculate the shortest achievable response time 
based on amplitude of the deviation and imposed con- 
strains of current state of the extended plant. 

[0081] These calculations can be made for complete envelope of 
imposed constraints. In example shown on Fig. 5 a lookup 
table for all specified nodes can be created. 

[0082] Complete block diagram of close-loop feedback operation 
is shown on Fig. 8. In response to the deviation driver 
functions are incremented byD. ft) and sensors readings 
are offset by S . (t). Implementation of controller algorithm 
may vary, and implement a set of optimal trajectories that 
satisfy different criteria (such as minimum peak power, 
fuel economy, fastest response, etc). 

[0083] Example of multithreaded MIMO controller 

[0084] This example illustrates construction of an apparatus that 
implements plurality of disclosed methods of this inven- 
tion. This example refers to the apparatus as a controller. 
It will become obvious after reading this description that 
various physical implementation of the controller are pos- 
sible due to wide availability and quick evolution of digital 
processing means and components for their creation. As 



an example controller with the same logical structure can 
be implemented using microcontroller, DSP, PDA, PC, and 
even discrete electronic components. Again it is obvious 
that some functional elements of the controller can be 
omitted at particular implementations to provide only 
minimal set of functionality required for specific solution. 
Again it is obvious that term multithreading can be ap- 
plied to a single digital processing device as well as to 
plurality of processing devices (commonly referred as 
multiprocessing), and although there is a difference in im- 
plementation details in these two cases the methods of 
the invention can be equally implemented. 
[0085] pig. 9 shows high level organization of the controller. Two 
real-time threads 501 and 502 are responsible for carry- 
ing out primary controller functions. Each of feed -forward 
or feedback pipe can be incorporated into control solution 
specific to particular application. Feed-forward pipe 501 
inputs receive control commands of desired plant outputs 
and produces plant control inputs. This pipe uses model 
data and current constraints. Some control commands 
may refer to use of specific time-optimal recipe. Note that 
time-optimal recipe could include various aspects of opti- 
mization, such as optimized for fuel efficiency, peak 



power consumption, etc. 501 requires minimum process- 
ing resources and depending on format of control com- 
mands and protocol can be implemented as a single exe- 
cution thread. It is more practical to separate command 
processing portion of 501 into separate thread to allow 
processing of more sophisticated control protocols and 
commands. 

[0086] Feedback pipe 502 inputs receives error signal that needs 
to be compensated. 502 uses recipes lookup 506 and 
scales down corresponding time-optimal control solution. 
Operations of feedback pipe 502 are very simple and re- 
quire minimal processing resources and it is practical to 
implement it in single real-time thread. 

[0087] Calibration pipe 508 operates in tie with calibration stack 
507. Pipe 508 receives set of control input sequences 
from model refinery 504. These sequences are executed 
in real time through output of the 508 pipe. Inputs of 508 
are acquired in real-time and pushed into FIFO stack 507 
along with time stamp. Because of simplicity of operation 
of 508 its both input and output operations may be per- 
formed by single real-time thread. 

[0088] | n some implementations it is possible to allocate single 
execution thread for 501 output and 502 and 508 opera- 



tions. 

[0089] Model refinery 504 implements an algorithm executing 
system identification method of the invention. It uses 
Model lookup table 505 to store models and retrieve 
model that require update. It is not required for 504 to 
operate in real-time. Due to complex nature of its opera- 
tions this module operates asynchronously and uses syn- 
chronization mechanisms when interacting with real-time 
threads. 

[0090] Time-optimal input factory 503 uses model data 505 and 
constraint lookup 509 to discover time-optimal sets of 
control solutions. These sets are stored in recipe lookup 
table 506. Factory 503 performs its operation in asyn- 
chronous mode and uses synchronization mechanisms to 
when interacting wit real-time threads. 

[° 091 ] State monitor 510 monitors current state of output and 
input of the plant and matches it to one of existing mod- 
els, which identifies current steady state of the plant. Op- 
erations of this model are time critical. Procedure of 
matching plant stimulus vector to matrix model can be 
quickly performed using SIMD processor instructions 
when vector has small dimensions (with today's mid range 
processors nearly real-time performance can be achieved 



for up to four dimensional vectors). Implementation algo- 
rithm for vectors of higher dimensions may restrict speed 
of module operations, but it will be sufficient for majority 
of applications that operates on sum megahertz rates. 

[0092] Possible implementations of the controller include: inte- 
gral controller that designed as single module containing 
all required elements; modular controller allows addition 
of new modules that expand its functions; and any inter- 
mediate implementation. 

[0093] Modular controller has an ability to operate with reduced 
set of features/operations some of them are illustrated 
through the following examples: calibration, close-loop 
feedback, open-loop feed-forward control, etc. 

[0094] Block diagram of calibration module of modular controller 
is shown on Fig 10. Referring to Fig 9 calibrator module 
includes 504, 505, 507, and 508 and executes the cali- 
bration methods described in previous sections of this in- 
vention. During operation it provides real-time control 
signals for each driver or/and channel controller module 
and records real-time sensor data. Particular form of real- 
time signals and data depends on specific of particular 
implementation details and may be realized as, but not 
limited, analog signals, serial digital input, parallel digital 



input or/and any type of data with time information that is 
sufficient for synchronization. Calibrator may be imple- 
mented as completely automated module that operates 
with no human intervention, or may be as generic as per- 
sonal computer with required set of peripheral devices. 
The essence of this invention does not depend on the 
form of practical implementation of this or any other 
module. The result of the calibrator operation a set of 
data described in previous sections. These data represent 
complete of incomplete models of the plant. The data 
recorded into long-term storage, which is required mod- 
ule for all features/operations of the controller. 

[0095] Long-term storage contains data sets of partitioned time 
and/or frequency series described in previous sections of 
the invention. Each set is sufficient for operation of a sin- 
gle channel of feedback 502 feed-forward 501 modules. 
Depending on practical realization of the controller these 
data sets may be stored in single long-term storage mod- 
ule or may be distributed on several smaller modules. 
Physical media for these long-term storage modules may 
be, but not limited, magnetic, optical, ROM, RAM, EPROM, 
FLASH and any other type of memory devices. 

[0096] input conditioner represents a set of similar and/or het- 



erogeneous devices that function as transducers convey- 
ing various real-world signals into electrical and/or digital 
form. An essential requirement for these transducers is 
fast conversion time. 

[0097] Mixer/ Equalizer converts all signals from transducers into 
digital or analog format. At the same time it may perform 
scaling of signal values to bring them to similar value 
range. This operation is not required if high precision 
arithmetic will be used in other digital modules, but it may 
help to increase computational precision. This module 
may mix different channels for various purposes 
(example: use redundant sensors to increase system relia- 
bility). Resulting number of output channels is always less 
or equal number of input channels, and exactly equal 
number of controller input channels. 

[0098] Driver modules convert digital and/or electrical signals 
into various forms of real-world energy. They provide 
stimulus that change current state of the plant. The num- 
ber of drivers is greater or equal to the number of control 
input channels of the plant, as example drivers may be 
redundant to increase reliability. 

[0099] close-loop feedback module is shown on Fig. 11. Its op- 
erations were described in previous sections. Referring to 



Fig. 9 this module encompasses operations of feedback 
pipe 502. 

[0100] Signals from Mixer/Equalizer in digital or analog format 
are transmitted to channel controller. Channel controller 
uses data from long-term storage module containing 
recipe lookup table 506. Output signal from this module 
controls operation of the driver or a set of redundant and/ 
or dependent drivers. The channel controller may use a 
user input that specifies a set of sensor reading or any 
other data that characterize the state of the plant. These 
data along with conditioned input data are used by state 
monitor 510. 

[0101] Block diagram of a open-loop feed-forward module is 

shown on Fig. 12. It operates according to description of 
feed-forward pipe 501 (refer to Fig. 9). The Automation 
control provides a set of command for channel controller 
that represent either optimal command stored in 506 
lookup table or custom input trajectory that based on 
model data 505 and current constraints obtained from 
509 lookup table. Channel controllers use data from long- 
term storage/s to determine output for driver in each 
control channel. The method that channel controllers use 
based on equation (eq.206). Creation of custom trajectory 



makes significant computational overhead for automation 
control when large number of channels is selected. 
[0102] For Automation control that uses optimal set of trajecto- 
ries all computations are done once by input factory 503 
(refer to Fig. 9), and results are stored in log-term stor- 
age. At run-time simple lookup and scaling of pre- 
computed D (t) is performed. 



